﻿#include <gtest/gtest.h>

#include "math/nonlinear_math.h"

#include <ctime>
#include <random>

using namespace CAGD;

#define PI 3.14159265358979323846

// 测试 diff
TEST(TestCalculus, diff)
{
    // Matrix operator()(FuncXX f, const Vector &x) const; // 多元向量函数的 Jacobi 矩阵
    // Matrix operator()(FuncX1 f, const Vector &x) const; // 多元函数的 Jacobi 矩阵
    // Vector operator()(Func1X f, double x) const;         // 一元向量函数的导向量
    // double operator()(Func11 f, double x) const;          // 一元函数的导数

    Func11 f1 = [](double x) { return sin(x); };
    Func1X f2 = [](double x) -> Vector { return {sin(x), cos(x)}; };
    FuncX1 f3 = [](Vector x) -> double { return 2 * sin(x[0]) + cos(x[1]); };
    FuncXX f4 = [](Vector x) -> Vector { return {2 * sin(x[0]) + cos(x[1]), sin(x[0]) - 2 * cos(x[1])}; };

    Jacobian J;
    auto r1 = J(f1, 0);
    auto r2 = J(f2, 0);
    auto r3 = J(f3, {0, 0});
    auto r4 = J(f4, {0, 0});

    EXPECT_DOUBLE_EQ(r1, 1);
    EXPECT_TRUE(r2 == Vector({1, 0}));
    EXPECT_TRUE(r3 == Matrix(1, 2, {2, 0}));
    EXPECT_TRUE(r4 == Matrix(2, 2, {2, 0, 1, 0}));
}

// 测试 newton-intergral
TEST(TestCalculus, newton_intergral)
{
    NewtonCotes I(8);
    Vector w = I.weights();

    auto f = [](double x) { return sin(PI * x); };
    EXPECT_FLOAT_EQ(I(f, 0, 1), 2 / PI);
}

// 测试 gauss-intergral
TEST(TestCalculus, gauss_intergral)
{
    GaussIntegrator I(8);
    Vector w = I.weights();
    Vector n = I.nodes();
    auto f = [](double x) { return sin(PI * x); };
    EXPECT_FLOAT_EQ(I(f, 0, 1), 2 / PI);
}

// 测试 NewtonSolver
TEST(TestNewtonSolver, solving)
{
    Func11 f1 = [](double x) {
        double y = x + sin(x) - 0.5 * x * x;
        return x;
    };

    NewtonSolver solver;
    double x0 = solver(f1, 0);
    EXPECT_DOUBLE_EQ(f1(x0), 0);
    Func11 f2 = [](double x) {
        double y = tan(x) + cos(x) - 0.5 * x * x;
        return x;
    };

    x0 = solver(f2, 0);
    EXPECT_DOUBLE_EQ(f2(x0), 0);

    FuncXX f3 = [=](Vector x) { return Vector{f1(x[0]), f2(x[1])}; };
    Vector X0(2);
    solver(f3, X0);
    EXPECT_TRUE(f3(X0) == Vector(2));

    FuncXX f4 = [=](Vector x) { return Vector{f1(f2(x[0])), f2(f1(x[1]))}; };
    Vector X1(2);
    solver(f4, X1);
    EXPECT_TRUE(f4(X1) == Vector(2));
}
